{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Optimal Power and Bandwidth Allocation in a Gaussian Channel\n",
    "by Robert Gowers, Roger Hill, Sami Al-Izzi, Timothy Pollington and Keith Briggs\n",
    "\n",
    "from Boyd and Vandenberghe, Convex Optimization, exercise 4.62 page 210\n",
    "\n",
    "Consider a system in which a central node transmits messages to $n$ receivers.  Each receiver channel $i \\in \\{1,...,n\\}$ has a transmit power $P_i$ and bandwidth $W_i$.  A fraction of the total power and bandwidth is allocated to each channel, such that $\\sum_{i=1}^{n}P_i = P_{tot}$ and $\\sum_{i=1}^{n}W_i = W_{tot}$.  Given some utility function of the bit rate of each channel, $u_i(R_i)$, the objective is to maximise the total utility $U = \\sum_{i=1}^{n}u_i(R_i)$.\n",
    "\n",
    "Assuming that each channel is corrupted by Gaussian white noise, the signal to noise ratio is given by $\\beta_i P_i/W_i$.  This means that the bit rate is given by:\n",
    "\n",
    "$R_i = \\alpha_i W_i \\log_2(1+\\beta_iP_i/W_i)$\n",
    "\n",
    "where $\\alpha_i$ and $\\beta_i$ are known positive constants.\n",
    "\n",
    "One of the simplest utility functions is the data rate itself, which also gives a convex objective function.\n",
    "\n",
    "The optimisation problem can be thus be formulated as:\n",
    "\n",
    "minimise $\\sum_{i=1}^{n}-\\alpha_i W_i \\log_2(1+\\beta_iP_i/W_i)$\n",
    "\n",
    "subject to $\\sum_{i=1}^{n}P_i = P_{tot} \\quad \\sum_{i=1}^{n}W_i = W_{tot} \\quad P \\succeq 0 \\quad W \\succeq 0$\n",
    "\n",
    "Although this is a convex optimisation problem, it must be rewritten in DCP form since $P_i$ and $W_i$ are variables and DCP prohibits dividing one variable by another directly.  In order to rewrite the problem in DCP format, we utilise the $\\texttt{kl_div}$ function in CVXPY, which calculates the Kullback-Leibler divergence.\n",
    "\n",
    "$\\text{kl_div}(x,y) = x\\log(x/y)-x+y$\n",
    "\n",
    "$-R_i = \\text{kl_div}(\\alpha_i W_i, \\alpha_i(W_i+\\beta_iP_i)) - \\alpha_i\\beta_iP_i$\n",
    "\n",
    "Now that the objective function is in DCP form, the problem can be solved using CVXPY."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#!/usr/bin/env python3\n",
    "# @author: R. Gowers, S. Al-Izzi, T. Pollington, R. Hill & K. Briggs\n",
    "\n",
    "import numpy as np\n",
    "import cvxpy as cvx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def optimal_power(n, a_val, b_val, P_tot=1.0, W_tot=1.0):\n",
    "  # Input parameters: α and β are constants from R_i equation\n",
    "  n=len(a_val)\n",
    "  if n!=len(b_val):\n",
    "    print('alpha and beta vectors must have same length!')\n",
    "    return 'failed',np.nan,np.nan,np.nan\n",
    "    \n",
    "  P=cvx.Variable(n)\n",
    "  W=cvx.Variable(n)\n",
    "  alpha=cvx.Parameter(n)\n",
    "  beta =cvx.Parameter(n)\n",
    "  alpha.value=np.array(a_val)\n",
    "  beta.value =np.array(b_val)\n",
    "  # This function will be used as the objective so must be DCP; \n",
    "  # i.e. elementwise multiplication must occur inside kl_div, not outside otherwise the solver does not know if it is DCP...\n",
    "  R=cvx.kl_div(cvx.mul_elemwise(alpha, W),\n",
    "               cvx.mul_elemwise(alpha, W + cvx.mul_elemwise(beta, P))) - \\\n",
    "    cvx.mul_elemwise(alpha, cvx.mul_elemwise(beta, P))\n",
    "  objective=cvx.Minimize(cvx.sum_entries(R))\n",
    "  constraints=[P>=0.0,\n",
    "               W>=0.0,\n",
    "               cvx.sum_entries(P)-P_tot==0.0,\n",
    "               cvx.sum_entries(W)-W_tot==0.0]\n",
    "  prob=cvx.Problem(objective, constraints)\n",
    "  prob.solve()\n",
    "  return prob.status,-prob.value,P.value,W.value"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Example\n",
    "\n",
    "Consider the case where there are 5 channels, $n=5$, $\\alpha = \\beta = (2.0,2.2,2.4,2.6,2.8)$, $P_{\\text{tot}} = 0.5$ and $W_{\\text{tot}}=1$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Status = optimal\n",
      "Optimal utility value = 2.451 \n",
      "Optimal power level:\n",
      " [[  1.150e-09]\n",
      " [  1.706e-09]\n",
      " [  2.754e-09]\n",
      " [  5.785e-09]\n",
      " [  5.000e-01]]\n",
      "Optimal bandwidth:\n",
      " [[  3.091e-09]\n",
      " [  3.956e-09]\n",
      " [  5.910e-09]\n",
      " [  1.193e-08]\n",
      " [  1.000e+00]]\n"
     ]
    }
   ],
   "source": [
    "  np.set_printoptions(precision=3)\n",
    "  n=5               # number of receivers in the system\n",
    "  a_val=np.arange(10,n+10)/(1.0*n)  # α\n",
    "  b_val=np.arange(10,n+10)/(1.0*n)  # β\n",
    "  P_tot=0.5\n",
    "  W_tot=1.0\n",
    "  status,utility,power,bandwidth=optimal_power(n,a_val,b_val,P_tot,W_tot)\n",
    "  print('Status: ',status)\n",
    "  print('Optimal utility value = %.4g '%utility)\n",
    "  print('Optimal power level:\\n', power)\n",
    "  print('Optimal bandwidth:\\n', bandwidth)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
